7-3 All Cause Mortality Day 28

Analyse of the all-cause mortality secondary outcome in ASCOT.

Authors
Affiliations
Rob Mahar

University of Melbourne

James Totterdell

University of Sydney

Published

July 13, 2022

Preamble

Load packages
library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(lubridate)
library(patchwork)
library(ggdist)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))
bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))
color_scheme_set("red")
options(digits = 4)
Prepare dataset
devtools::load_all()
all_dat <- read_all_no_daily()

# Anticoagulation set (intention-to-treat, including missing)
acs_itt_dat <- all_dat %>% 
  filter_acs_itt() %>%
  transmute_model_cols_grp_aus_nz()

# Anticoagulation set (intention-to-treat, excluding missing)
acs_itt_nona_dat <- acs_itt_dat %>% 
  filter(!is.na(D28_death))  

# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
  filter_concurrent_intermediate()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>% 
  filter(!is.na(D28_death))

# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
  filter_concurrent_std_aspirin()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>% 
  filter(!is.na(D28_death))

# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
  filter_concurrent_therapeutic()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>% 
  filter(!is.na(D28_death))
Load the Stan models
model_full <- cmdstan_model(file.path("stan", "binary", "logistic_site_epoch.stan"))
model_site_only <- cmdstan_model(file.path("stan", "binary", "logistic_site.stan"))
model_epoch_only <- cmdstan_model(file.path("stan", "binary", "logistic_epoch.stan")) 
model_fixed_only <- cmdstan_model(file.path("stan", "binary", "logistic.stan"))
Helper functions for this analysis script
make_stan_data <- function(
    dat, 
    vars = NULL, 
    outcome = NULL, 
    beta_sd = NULL, 
    ctr = contr.orthonorm) {

  X <- make_X_design(dat, vars, ctr)
  nXassign <- sum(grepl("rand", colnames(X))) - 1
  beta_sd <- c(2.5, rep(1, nXassign), beta_sd)

  y <- dat[[outcome]]
  N <- dim(X)[1]
  K <- dim(X)[2]  
  
    epoch  <- dat$epoch
  M_epoch  <- max(dat$epoch)
    region <- dat[["ctry_num"]]
  M_region <- max(region)
    site <- dat[["site_num"]]
  M_site <- max(site)
  region_by_site <- region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)
  
  list(N = N, K = K, X = X, y = y,
       M_region = M_region, region = region,
       M_site = M_site, site = site,
       M_epoch = M_epoch, epoch = epoch,
       region_by_site = region_by_site,
       beta_sd = beta_sd)
}

make_stan_data_concurrent <- function(
    dat, 
    vars = NULL, 
    outcome = NULL, 
    beta_sd = NULL, 
    ctr = contr.orthonorm) {
  
  X <- make_X_design(dat, vars, ctr, includeA = FALSE)
  nXassign <- sum(grepl("rand", colnames(X))) - 1
  beta_sd <- c(2.5, rep(1, nXassign), beta_sd)
  y <- dat[[outcome]]
  N <- dim(X)[1]
  K <- dim(X)[2]  
  list(N = N, K = K, X = X, y = y, beta_sd = beta_sd)
}

make_mortality_table <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3", "C4"),
      labels = str_replace(intervention_labels()$CAssignment[-1], "<br>", " "))) %>%
    group_by(CAssignment) %>%
    summarise(
      Patients = n(),
      Known = sprintf(
        "%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
      Missing = sprintf(
        "%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
      `Died by day 28` = sprintf(
        "%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Patients = n(),
          Known = sprintf(
            "%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
          Missing = sprintf(
            "%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
          `Died by day 28` = sprintf(
            "%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(key, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(font_size = 10, latex_options = "HOLD_position")  %>%
    row_spec(nrow(tab), bold = TRUE)
}

Descriptive

Anticoagulation

Summary table
save_tex_table(make_mortality_table(
      all_dat %>%
      filter_acs_itt(),
      "latex"), 
      "outcomes/secondary/anticoagulation-summary-mortality")
make_mortality_table(
    all_dat %>%
      filter_acs_itt()
)
n (\%) Patients Known Missing Died by day 28
Low dose 610 596 (97.7) 14 (2.3) 19 (3.2)
Intermediate dose 613 603 (98.4) 10 (1.6) 15 (2.5)
Low dose with aspirin 283 281 (99.3) 2 (0.7) 10 (3.6)
Therapeutic dose 50 50 (100.0) 0 (0.0) 6 (12.0)
Overall 1556 1530 (98.3) 26 (1.7) 50 (3.3)
Summary of mortality to day 28 by anticoagulation intervention.

Age

Relationship age to outcome
agedat <- acs_itt_dat %>%
  dplyr::count(D28_death, AgeAtEntry) %>% 
  spread(D28_death, n, fill = 0) %>% 
  mutate(
    n = `0` + `1` + `<NA>`,
    p = `1` / (`1` + `0`))
agemod <- glm(
  cbind(`1`, `0`) ~ AgeAtEntry, 
  data = agedat, 
  family = binomial())
agedat <- agedat %>%
  mutate(
    ypred = predict(agemod, newdata = agedat, type = "response")
  )
p1 <- ggplot(agedat, aes(AgeAtEntry, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry")
p2 <- ggplot(agedat, aes(AgeAtEntry, p)) +
    geom_point() +
    geom_vline(xintercept = 60, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\ndied by day 28", x = "Age at entry")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-age.pdf"), p, height = 2, width = 6)
p

Country

Relationship country to outcome
tdat <- all_dat %>%
  filter_acs_itt() %>%
  dplyr::count(Country = fct_infreq(PT_CountryName), D28_death) %>%
  group_by(Country) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(tdat, aes(Country, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Country of enrolment")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-country.pdf"), p, height = 2, width = 6)
p

Site

Relationship site to outcome
tdat <- all_dat %>%
  filter_acs_itt() %>%
  dplyr::count(
    Country_lab = Country,
    Site_lab = fct_infreq(Location),
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")),
    Site = PT_LocationName,
    D28_death) %>%
  group_by(Country, Site) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  ungroup()

p1 <- ggplot(tdat, aes(Site_lab, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(tdat, aes(Site_lab, p_1)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Site of enrolment") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p <- p1 / p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-country-site.pdf"), p, height = 4, width = 6.25)
p

Calendar Time

Relationship calendar date to outcome
caldat <- all_dat %>% 
  filter_acs_itt() %>%
  dplyr::count(D28_death, yr = year(RandDate), mth = month(RandDate)) %>% 
  spread(D28_death, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`),
         n = `1` + `0` + `<NA>`)
p1 <- ggplot(caldat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12) +
  ylim(0, 0.15)
p2 <- ggplot(caldat, aes(mth, n)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p <- p2 | p1
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-calendar-time.pdf"), p, height = 2, width = 6)
p

Modelling

The analysis for all cause mortality is specified equivalently to that for the primary outcome. It includes intervention as randomised, age category, country, site nested within country, epoch, and intervention eligibility. The main analysis restricts the population to only those participants randomised to the anticoagulation domain.

Primary model

Data and sampling
seed <- 59579
mdat <- make_stan_data(
    dat = acs_itt_nona_dat, 
    vars    = c("inelgc3", "agegte60", "ctry"),
    outcome = "D28_death",
    beta_sd = c(10, 2.5, 1, 1))
snk <- capture.output(
  mfit <- model_full$sample(
    data = mdat,
    seed = seed,
    refresh = 0,
    iter_warmup = 1000,
    iter_sampling = 2500,
    chains = 8,
    parallel_chains = 8,
    adapt_delta = 0.99
))
mfit$save_object(file.path("outputs", "models", "secondary", "7-3.rds"))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(mdat$X)
names(mdrws$gamma_epoch) <- acs_itt_nona_dat %>% dplyr::count(epoch, epoch_lab) %>% pull(epoch_lab)
names(mdrws$gamma_site) <- acs_itt_nona_dat %>% dplyr::count(site_num, site) %>% pull(site)
site_facet <- factor(mdat$region_by_site, labels = c("India", "Australia\nNew Zealand", "Nepal"))

mdrws$Acon <- attr(mdat$X, "contrasts")$randA %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))]
mdrws$Atrt <- mdrws$Acon[-1] - mdrws$Acon[1]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c(
  "Intermediate vs low" = exp(mdrws$Ctrt[1]),
  "Low with aspirin vs low" = exp(mdrws$Ctrt[2]),
  "Therapeutic vs low" = exp(mdrws$Ctrt[3]),
  "Intermediate vs low with aspirin" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]),
  "Intermediate vs therapeutic" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[3]),
  "Low with aspirin vs therapeutic" = exp(mdrws$Ctrt[2] - mdrws$Ctrt[3])
)
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic")),
  setNames(exp(mdrws$beta[7:8]), c("Ineligible aspirin", "Age \u2265 60")),
  setNames(exp(mdrws$beta[9:10]), c("Australia/New Zealand", "Nepal"))
)
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-3-primary-model-acs-itt-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.84 (0.42, 1.70) 0.89 (0.33) 0.69
Low with aspirin 0.82 (0.35, 1.84) 0.89 (0.39) 0.69
Therapeutic 3.17 (0.96, 10.05) 3.77 (2.43) 0.03
Ineligible aspirin 4.73 (1.23, 16.96) 5.82 (4.21) 0.01
Age ≥ 60 2.08 (1.12, 3.86) 2.19 (0.72) 0.01
Australia/New Zealand 0.62 (0.13, 2.75) 0.82 (0.73) 0.73
Nepal 2.85 (0.57, 11.08) 3.61 (2.87) 0.09
Epoch and site terms
p <- plot_epoch_site_terms(
  exp(mdrws$gamma_epoch), 
  exp(mdrws$gamma_site), 
  site_facet
)
pth <- file.path(
  "outputs", "figures", "outcomes", "secondary", "7-3-primary-model-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Code
p <- plot_or_densities(mdrws$compare)
p

Code
mfit$summary(variables = c("beta", "gamma_site", "gamma_epoch"))
# A tibble: 48 × 10
   variable    mean  median    sd   mad     q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>  <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -4.73   -4.68   0.765 0.754 -6.05  -3.56    1.00   11020.   13326.
 2 beta[2]   0.975   0.977  0.483 0.482  0.189  1.78    1.00   24918.   15667.
 3 beta[3]  -0.222  -0.220  0.282 0.281 -0.689  0.239   1.00   26301.   16533.
 4 beta[4]  -0.499  -0.501  0.330 0.326 -1.04   0.0440  1.00   24973.   16489.
 5 beta[5]  -0.999  -1.00   0.828 0.822 -2.37   0.364   1.00   26126.   15510.
 6 beta[6]  -0.0637 -0.0680 0.845 0.854 -1.43   1.31    1.00   31265.   15003.
 7 beta[7]   1.54    1.55   0.664 0.663  0.449  2.61    1.00   27291.   15723.
 8 beta[8]   0.733   0.734  0.318 0.318  0.211  1.26    1.00   28048.   15306.
 9 beta[9]  -0.498  -0.485  0.792 0.802 -1.81   0.786   1.00   21930.   15612.
10 beta[10]  1.02    1.05   0.748 0.719 -0.265  2.20    1.00   15008.   13871.
# … with 38 more rows
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7989 0.8661 0.8909 0.8373 0.8476 0.8309 0.8887 0.8998
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["gamma_site"])

Code
mcmc_trace(mdrws["gamma_epoch"])

Code
mcmc_trace(mdrws[c("tau_site", "tau_epoch")])

Posterior predictive

Code
y_ppc <- mdrws$y_ppc
ppc_dat <- bind_cols(acs_itt_nona_dat, tibble(y_ppc = y_ppc))

grp_ppc <- function(grp) {
  ppc_dat %>%
    group_by(grp = !!grp) %>%
    summarise(
      mean_y = mean(D28_death),
      rvar_mean_y_ppc = rvar_mean(y_ppc)
    )
}
plot_grp_ppc <- function(dat, lab = "") {
  dat %>%
    ggplot(aes(y = grp, xdist = rvar_mean_y_ppc)) +
    stat_interval(size = 2) +
    geom_point(aes(x = mean_y, y = 1:nrow(dat)), colour = "red", shape = 23) +
    labs(
      x = "Posterior predictive\nproportion", 
      y = lab)  
}

ppc_C <- grp_ppc(sym("CAssignment"))
ppc_ctry <- grp_ppc(sym("country"))
ppc_epoch <- grp_ppc(sym("epoch"))
ppc_site <- ppc_dat %>%
  group_by(Country = country, grp = site_raw) %>%
  summarise(
    mean_y = mean(D28_death),
    rvar_mean_y_ppc = rvar_mean(y_ppc)
  )

p1 <- plot_grp_ppc(ppc_C, "Anticoagulation\nintervention") + labs(x = "")
p2 <- plot_grp_ppc(ppc_ctry, "Country") + labs(x = "")
p3 <- plot_grp_ppc(ppc_epoch, "Epoch") + labs(x = "")
p4 <- plot_grp_ppc(ppc_site %>% filter(Country == "IN"), "Sites\nIndia")  + labs(x = "")
p5 <- plot_grp_ppc(ppc_site %>% filter(Country == "AU"), "Sites\nAustralia") + labs(x = "")
p6 <- plot_grp_ppc(ppc_site %>% filter(Country == "NP"), "Sites\nNepal")
p7 <- plot_grp_ppc(ppc_site %>% filter(Country == "NZ"), "Sites\nNZ")
p <- ((p3 | p1 / p2) / 
        ( (p4 | p5) / (p6 | p7) + 
            plot_layout(heights = c(3, 1)) ) ) +
  plot_layout(heights = c(1, 1.5), guides = "collect")
pth <- file.path("outputs", "figures", "outcomes", "secondary",
                 "7-3-primary-model-acs-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.5)
p

Sensitivity: Reduced model (no site terms)

The primary model is re-fit without the hierarchical site terms.

Fit model
mdat <- make_stan_data(
  acs_itt_nona_dat,
  c("inelgc3", "agegte60", "ctry"),
  "D28_death",
  c(10, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_epoch_only$sample(
    data = mdat,
    seed = 14646,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0,
    adapt_delta = 0.99
))
mdrws <- as_draws_rvars(mfit$draws(c("beta", "gamma_epoch", "tau_epoch")))
names(mdrws$beta) <- colnames(mdat$X)
epoch_map <- acs_itt_nona_dat %>% dplyr::count(epoch, epoch_lab)
names(mdrws$gamma_epoch) <- epoch_map$epoch_lab
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:4]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic")),
  setNames(exp(mdrws$beta[-(1:6)]), c("Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal"))
)
Odds ratio summary
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.83 (0.42, 1.61) 0.88 (0.31) 0.71
Low with aspirin 0.84 (0.38, 1.83) 0.91 (0.37) 0.67
Therapeutic 3.39 (1.04, 10.68) 4.04 (2.61) 0.02
Ineligible aspirin 2.86 (0.85, 8.26) 3.30 (1.97) 0.05
Age ≥ 60 2.32 (1.28, 4.25) 2.44 (0.76) 0.00
Australia/New Zealand 0.76 (0.16, 3.05) 0.97 (0.78) 0.65
Nepal 4.56 (1.55, 13.75) 5.34 (3.29) 0.00
Odds ratio summary table for no epoch model.
Code
plot_epoch_terms(exp(mdrws$gamma_epoch))

Code
mfit$summary(variables = c("beta", "gamma_epoch", "tau_epoch"))
# A tibble: 23 × 10
   variable    mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -5.01   -4.97   0.663 0.652 -6.16   -3.99    1.00   14715.   12970.
 2 beta[2]   1.01    1.00   0.478 0.478  0.224   1.78    1.00   20283.   15713.
 3 beta[3]  -0.246  -0.246  0.274 0.274 -0.701   0.201   1.00   22550.   15654.
 4 beta[4]  -0.548  -0.547  0.326 0.330 -1.08   -0.0154  1.00   19992.   15215.
 5 beta[5]  -0.936  -0.932  0.811 0.803 -2.30    0.386   1.00   25416.   16103.
 6 beta[6]  -0.0529 -0.0600 0.856 0.852 -1.45    1.37    1.00   26701.   13969.
 7 beta[7]   1.03    1.05   0.585 0.581  0.0344  1.96    1.00   25791.   13258.
 8 beta[8]   0.844   0.843  0.303 0.303  0.349   1.34    1.00   27646.   15597.
 9 beta[9]  -0.306  -0.278  0.756 0.753 -1.59    0.900   1.00   24776.   15176.
10 beta[10]  1.52    1.52   0.557 0.553  0.602   2.43    1.00   15983.   14684.
# … with 13 more rows
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8841 0.9230 0.8853 0.8829 0.9267 0.8939 0.8181 0.8834
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["gamma_epoch"])

Code
mcmc_trace(mdrws["tau_epoch"])

Sensitivity: Reduced model (no epoch terms)

The primary model is re-fit without the epoch terms.

Fit model
mdat <- make_stan_data(
  acs_itt_nona_dat,
  c("inelgc3", "agegte60", "ctry"),
  "D28_death",
  c(10, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_site_only$sample(
    data = mdat,
    seed = 123513,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0,
    adapt_delta = 0.99
))

mdrws <- as_draws_rvars(mfit$draws(c("beta", "gamma_site", "tau_site")))
site_map <- acs_itt_nona_dat %>% dplyr::count(site_num, site)
site_facet <-   factor(mdat$region_by_site, labels = c("India", "Australia\nNew Zealand", "Nepal"))
names(mdrws$beta) <- colnames(mdat$X)
names(mdrws$gamma_site) <- site_map$site
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:4]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic")),
  setNames(exp(mdrws$beta[-(1:6)]), c("Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal"))
)
Odds ratio summary
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.84 (0.42, 1.67) 0.89 (0.33) 0.69
Low with aspirin 0.88 (0.39, 1.95) 0.96 (0.41) 0.62
Therapeutic 2.93 (0.90, 9.11) 3.46 (2.22) 0.04
Ineligible aspirin 5.42 (1.37, 19.17) 6.63 (4.78) 0.01
Age ≥ 60 1.97 (1.07, 3.65) 2.07 (0.67) 0.01
Australia/New Zealand 0.59 (0.12, 2.64) 0.80 (0.71) 0.75
Nepal 2.28 (0.57, 7.75) 2.78 (2.03) 0.11
Odds ratio summary table for no epoch model.
Code
plot_site_terms(exp(mdrws$gamma_site), site_facet)

Code
mfit$summary(variables = c("beta", "gamma_site", "tau_site"))
# A tibble: 39 × 10
   variable    mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -4.07   -4.05   0.519 0.504 -4.95   -3.26    1.00    8867.   12171.
 2 beta[2]   0.859   0.858  0.471 0.467  0.0896  1.63    1.00   22633.   14901.
 3 beta[3]  -0.221  -0.218  0.281 0.278 -0.685   0.240   1.00   26235.   16195.
 4 beta[4]  -0.505  -0.503  0.326 0.322 -1.05    0.0311  1.00   24775.   15954.
 5 beta[5]  -1.04   -1.04   0.834 0.841 -2.41    0.323   1.00   25838.   16057.
 6 beta[6]  -0.0811 -0.0795 0.849 0.844 -1.47    1.30    1.00   30912.   13477.
 7 beta[7]   1.67    1.69   0.667 0.654  0.568   2.75    1.00   23886.   15322.
 8 beta[8]   0.679   0.679  0.313 0.314  0.162   1.20    1.00   29941.   14381.
 9 beta[9]  -0.531  -0.521  0.790 0.786 -1.85    0.752   1.00   22361.   15767.
10 beta[10]  0.807   0.826  0.666 0.650 -0.323   1.86    1.00   13402.   14508.
# … with 29 more rows
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9495 0.8738 0.9165 0.8951 0.9422 0.8839 0.9512 0.8555
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["gamma_site"])

Code
mcmc_trace(mdrws["tau_site"])

Sensitivity: Reduced model (no site or epoch)

The primary model is re-fit without the site and epoch terms.

Fit reduced model (no site or epoch)
mdat <- make_stan_data(
  acs_itt_nona_dat,
  c("inelgc3", "agegte60", "ctry"),
  "D28_death",
  c(10, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 12415,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0,
    adapt_delta = 0.99
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:4]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic")),
  setNames(exp(mdrws$beta[-(1:6)]), c("Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal"))
)
Odds ratio summary
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.81 (0.41, 1.58) 0.86 (0.30) 0.73
Low with aspirin 1.26 (0.57, 2.68) 1.35 (0.55) 0.28
Therapeutic 2.58 (0.84, 7.33) 2.96 (1.73) 0.05
Ineligible aspirin 4.20 (1.28, 11.39) 4.75 (2.67) 0.01
Age ≥ 60 2.36 (1.30, 4.19) 2.46 (0.76) 0.00
Australia/New Zealand 0.57 (0.13, 2.18) 0.72 (0.55) 0.79
Nepal 1.97 (0.81, 4.39) 2.13 (0.93) 0.06
Odds ratio summary table for no epoch model.
Code
mfit$summary(variables = c("beta"))
# A tibble: 10 × 10
   variable    mean  median    sd   mad      q5    q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -3.65   -3.65   0.257 0.252 -4.08   -3.23   1.00   11510.   13122.
 2 beta[2]   0.522   0.524  0.443 0.442 -0.215   1.24   1.00   13728.   15091.
 3 beta[3]  -0.274  -0.273  0.267 0.264 -0.711   0.167  1.00   16229.   13075.
 4 beta[4]  -0.631  -0.631  0.313 0.312 -1.14   -0.115  1.00   15687.   14173.
 5 beta[5]  -1.12   -1.11   0.794 0.783 -2.44    0.188  1.00   17669.   15548.
 6 beta[6]  -0.0715 -0.0679 0.831 0.815 -1.44    1.28   1.00   21121.   14415.
 7 beta[7]   1.41    1.43   0.559 0.552  0.451   2.28   1.00   18010.   14037.
 8 beta[8]   0.855   0.858  0.301 0.303  0.360   1.34   1.00   16662.   14081.
 9 beta[9]  -0.580  -0.566  0.725 0.729 -1.80    0.576  1.00   18319.   14119.
10 beta[10]  0.665   0.680  0.432 0.429 -0.0621  1.36   1.00   14109.   14057.
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9655 1.0383 0.9775 1.0953 1.0490 0.9577 1.0132 1.0193
Code
mcmc_trace(mdrws["beta"])

Sensitivity: Concurrent Controls

A sensitivity analysis is performed on the reduced concurrent control analysis sets for each intervention.

Intermediate

  • Set: ACS-ITT-intermediate
  • Covariates: anticoagulation intervention, age group, region of enrolment.
Overview of outcome
save_tex_table(
  make_mortality_table(acs_itt_concurc2_dat, "latex"),
  "outcomes/secondary/7-3-anticoagulation-concurrent-intermediate-summary")
make_mortality_table(acs_itt_concurc2_dat)
n (\%) Patients Known Missing Died by day 28
Low dose 610 596 (97.7) 14 (2.3) 19 (3.2)
Intermediate dose 613 603 (98.4) 10 (1.6) 15 (2.5)
Overall 1223 1199 (98.0) 24 (2.0) 34 (2.8)
Summary of outcome for intermediate dose concurrent enrolments.
Fit model
mdat <- make_stan_data_concurrent(
  acs_itt_concurc2_nona_dat,
  c("agegte60", "ctry"),
  "D28_death",
  c(2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 38135,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0,
    adapt_delta = 0.95
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]))
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate")),
  setNames(exp(mdrws$beta[-(1:2)]), c("Age \u2265 60", "Australia/New Zealand", "Nepal"))
)
Odds ratio summary
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-3-primary-model-acs-itt-concurrent-intermediate-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.78 (0.39, 1.54) 0.83 (0.29) 0.76
Age ≥ 60 2.84 (1.44, 5.62) 3.01 (1.08) 0.00
Australia/New Zealand 0.46 (0.12, 1.41) 0.54 (0.34) 0.91
Nepal 1.31 (0.46, 3.22) 1.45 (0.72) 0.30
Odds ratio summary table for model.
Posterior contrast
p <- plot_or_densities(mdrws$compare)
p

Code
mfit$summary(variables = c("beta"))
# A tibble: 5 × 10
  variable   mean median    sd   mad     q5    q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -3.95  -3.94  0.258 0.256 -4.39  -3.54   1.00   11239.   11056.
2 beta[2]  -0.175 -0.173 0.246 0.246 -0.582  0.226  1.00   18118.   13824.
3 beta[3]   1.04   1.04  0.346 0.344  0.473  1.61   1.00   12024.   11313.
4 beta[4]  -0.804 -0.766 0.633 0.625 -1.89   0.169  1.00   17707.   12582.
5 beta[5]   0.250  0.269 0.497 0.490 -0.601  1.04   1.00   16895.   13100.
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9724 1.0293 1.0048 0.9900 1.0677 1.0275 1.0277 1.0140
Code
mcmc_trace(mdrws["beta"])

Low-dose with aspirin

  • Set: ACS-ITT-aspirin
  • Covariates: anticoagulation intervention, age group, region of enrolment.
Overview of outcome
save_tex_table(
  make_mortality_table(acs_itt_concurc3_dat, "latex"),
  "outcomes/secondary/7-3-anticoagulation-concurrent-stdaspirin-summary")
make_mortality_table(acs_itt_concurc3_dat)
n (\%) Patients Known Missing Died by day 28
Low dose 299 291 (97.3) 8 (2.7) 11 (3.8)
Intermediate dose 298 293 (98.3) 5 (1.7) 12 (4.1)
Low dose with aspirin 283 281 (99.3) 2 (0.7) 10 (3.6)
Overall 880 865 (98.3) 15 (1.7) 33 (3.8)
Summary of outcome for low-dose with aspirin concurrent enrolments.
Fit model
mdat <- make_stan_data_concurrent(
  acs_itt_concurc3_nona_dat,
  c("agegte60", "ctry"),
  "D28_death",
  c(2.5, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 13578,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:3]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Low with aspirin vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs low with aspirin" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin")),
  setNames(exp(mdrws$beta[-(1:3)]), c("Age \u2265 60", "Australia/New Zealand"))
)
Odds ratio summary
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-3-primary-model-acs-itt-concurrent-stdaspirin-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 1.09 (0.48, 2.48) 1.19 (0.53) 0.42
Low with aspirin 0.89 (0.38, 2.07) 0.97 (0.44) 0.61
Age ≥ 60 2.45 (1.22, 4.87) 2.60 (0.95) 0.01
Australia/New Zealand 0.41 (0.08, 1.58) 0.52 (0.40) 0.89
Odds ratio summary table for model.
Posterior contrast
p <- plot_or_densities(mdrws$compare)
p

Code
mfit$summary(variables = c("beta"))
# A tibble: 5 × 10
  variable    mean  median    sd   mad     q5    q95  rhat ess_bulk ess_tail
  <chr>      <dbl>   <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -3.57   -3.56   0.246 0.244 -3.99  -3.18   1.00   16232.   14278.
2 beta[2]  -0.148  -0.145  0.301 0.301 -0.649  0.338  1.00   22800.   14161.
3 beta[3]   0.0127  0.0187 0.300 0.299 -0.493  0.495  1.00   22693.   13702.
4 beta[4]   0.894   0.895  0.355 0.355  0.311  1.47   1.00   16812.   14811.
5 beta[5]  -0.913  -0.881  0.752 0.746 -2.21   0.269  1.00   22255.   13232.
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.097 1.026 1.149 1.007 1.063 1.109 1.036 1.038
Code
mcmc_trace(mdrws["beta"])

Therapeutic

  • Set: ACS-ITT-therapeutic
  • Covariates: anticoagulation intervention, age group, region of enrolment.
Overview of outcome
save_tex_table(
  make_mortality_table(acs_itt_concurc4_dat, "latex"),
  "outcomes/secondary/7-3-anticoagulation-concurrent-therapeutic-summary")
make_mortality_table(acs_itt_concurc4_dat)
n (\%) Patients Known Missing Died by day 28
Low dose 79 75 (94.9) 4 (5.1) 3 (4.0)
Intermediate dose 65 63 (96.9) 2 (3.1) 1 (1.6)
Therapeutic dose 50 50 (100.0) 0 (0.0) 6 (12.0)
Overall 194 188 (96.9) 6 (3.1) 10 (5.3)
Summary of outcome for therapeutic dose concurrent enrolments.
Fit model
mdat <- make_stan_data_concurrent(
  acs_itt_concurc4_nona_dat,
  c("agegte60", "ctry"),
  "D28_death",
  c(2.5, 1, 1)
)
snk <- capture.output(
  mfit <- model_fixed_only$sample(
    data = mdat,
    seed = 13578,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws(c("beta")))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[2:3]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Therapeutic vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs therapeutic" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c(
  setNames(exp(mdrws$Ctrt), c("Intermediate", "Therapeutic")),
  setNames(exp(mdrws$beta[-(1:3)]), c("Age \u2265 60", "India", "Australia/New Zealand"))
)
Odds ratio summary
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-3-primary-model-acs-itt-concurrent-therapeutic-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.66 (0.13, 3.11) 0.91 (0.86) 0.69
Therapeutic 2.64 (0.70, 10.68) 3.38 (2.75) 0.07
Age ≥ 60 2.80 (0.78, 11.41) 3.60 (3.05) 0.06
India 0.42 (0.08, 1.81) 0.55 (0.46) 0.87
Australia/New Zealand 0.28 (0.06, 1.02) 0.35 (0.26) 0.97
Odds ratio summary table for model.
Posterior contrast
p <- plot_or_densities(mdrws$compare)
p

Code
mfit$summary(variables = c("beta"))
# A tibble: 6 × 10
  variable   mean median    sd   mad      q5    q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl> <dbl> <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -3.33  -3.29  0.566 0.551 -4.32   -2.47   1.00   12097.   10881.
2 beta[2]   0.986  0.969 0.555 0.543  0.0946  1.93   1.00   18013.   14124.
3 beta[3]  -0.228 -0.220 0.526 0.520 -1.10    0.622  1.00   19625.   14656.
4 beta[4]   1.04   1.03  0.676 0.660 -0.0390  2.17   1.00   12790.   10806.
5 beta[5]  -0.885 -0.861 0.784 0.779 -2.22    0.362  1.00   20168.   13677.
6 beta[6]  -1.29  -1.26  0.711 0.711 -2.50   -0.179  1.00   20226.   13969.
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0004 1.0536 1.0274 1.0999 1.0111 0.9799 1.0311 0.9898
Code
mcmc_trace(mdrws["beta"])